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Study of near-surface models in large-eddy 
simulations of a neutrally stratified atmospheric 

boundary layer 
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1. Motivation and objectives 

Large-eddy simulation (LES) is a widely used technique in atmospheric modeling re- 
search, partly because of the difficulties involved in observational studies and field exper- 
iments to obtain information about the turbulent structure of the atmosphere (Stevens 
& Lenschow 2001). In LES, large, unsteady, three dimensional structures are resolved 
and small structures that are not resolved on the computational grid are modeled. A 
filtering operation is applied to distinguish between resolved and unresolved scales. Un- 
resolved motions are believed to be universal, and simple models should be sufficient to 
parameterize them, provided that a major fraction of the energetic large scales are re- 
solved by the spatial and temporal resolution (Sagaut 2002). However, this requirement 
becomes stringent to fulfill in the close proximity of the surface/wall. Traditionally, the 
no-slip boundary has been referred to as “wall” in engineering context, and as “surface” 
in atmospheric science context. We will use the latter terminology throughout the paper. 

As the surface is approached anisotropy in turbulence structure increases and the 
length scale of the flow structures diminish rapidly, requiring too fine a spatial and a 
temporal resolution to numerically resolve a large fraction of the energetic scales. Hence, 
fully resolved LES of atmospheric or any other wall-bounded high Reynolds number 
flow would be extremely expensive in terms of computational resources. In fact, it is 
this obstacle that motivates surface/wall modeling research in LES. Furthermore, the 
roughness of the surface underlying the atmospheric boundary layer has always been a 
complicated issue, and needs to be considered in modeling too. 

In surface modeling, no-slip boundary conditions are not applied directly because the 
implied stress would be overestimated on a coarse grid. The simplest surface model 
that has long been adopted in atmospheric boundary layer simulations assumes that the 
logarithmic law holds within the surface layer, and stresses are imposed as boundary con- 
ditions at the surface(Schumann 1975; Moeng 1986). Cabot, Jimenez & Baggett (1999) 
showed that such models have problems for high Reynolds number LES simulations with 
coarse numerical resolution. Cabot (1997) indicated that providing accurate mean sur- 
face stresses is not sufficient to overcome the overall poor predictions. He suggested that, 
given a coarse numerical resolution, typical subgrid-scale (SGS) models used in LES do 
not perform well to predict Reynolds stresses in the near-surface region accurately. 

Boundary layer approximations have been proposed to extend the application of sur- 
face models for separating flows with adverse pressure gradient, (Balaras et al. 1996; 
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Cabot & Moin 2000). In this approach, a simplified set of turbulent boundary layer 
equations, adopting a Reynolds-averaged Navier-Stokes (RANS) type eddy viscosity, are 
solved on an embedded mesh near the wall. The surface stress is then calculated from the 
computed velocity profile. Within this framework, Wang & Moin (1992) utilized a dy- 
namically adjusted mixing length eddy viscosity, and showed that their model performs 
significantly better than the simpler surface modeling approaches. Surface modeling for 
LES continues to be an active area of research for high Reynolds number flow, coarse 
resolution simulations. Latest research in this area focuses on incorporating suboptimal 
control theory to improve the predictive capability of the surface models (Templeton, 
Wang & Moin 2002). 

In atmospheric boundary layer simulations, the practice has been to employ LES mod- 
els away from the surface and make a transition towards ensemble-averaged (RANS) 
models as the surface is approached. This idea has been adopted through different for- 
mulations. For instance, Sullivan et al. (1994) proposed an eddy viscosity model in which 
the so-called isotropy factor controls the transition from LES to a RANS type simu- 
lation, and accounts for the anisotropy effects at the same time. Mason & Thomson 
(1992) introduced a modified length scale that is a matching between the LES filter size 
and the distance from the surface. Brown, Hobson & Wood (2001) suggested using a 
canopy model with the aim of overcoming the difficulties in surface modeling. Tradition- 
ally, canopy models have been adopted to study the effect of flow within the canopy 
layer. In this approach, instead of modifying the eddy viscosity or the length scale, ad- 
ditional turbulent stresses are added to the turbulent stresses that are already modeled 
by the LES SGS model. Chow & Street (2002) and Kirkpatrick et al. (2003) adopted 
the canopy model of Brown et al. (2001) for atmospheric boundary layer simulations, 
and reported improved predictions. We will refer to these models, in which additional 
stresses are added onto sugrid-scale models to improve the predictions in the vicinity of 
the surface, as “near-surface” models. We reserve the term “surface model” to include 
both the models for surface boundary condition and the near-surface models. 

In what follows, we present two near-surface models that have found use in atmp- 
spheric modeling. We also suggest a simpler eddy viscosity model that adopts Prandtl’s 
mixing length model (Prandtl 1925) in the vicinity of the surface and blends in with the 
dynamic Smagorinsky model (Germano et al. 1991) away from the surface. We evaluate 
the performance of these surface model by simulating a neutrally stratified atmospheric 
boundary layer. 


2. Governing equations 

The governing equations for LES of a neutrally stratified atmospheric boundary layer 
are the filtered Navier-Stokes equations 



duj d(ujUj) _ dp_ _ drjj 

dt dxj dxi jUk dxj 


(2.1) 


where turbulent stresses are defined as nj = uiuj — UiUj , fj is the Coriolis parameter, 
eijk is the permutation tensor, u and p are the filtered velocity and filtered dynamic 
pressure, respectively. 

The lateral boundary conditions are periodic and a stress- free condition is imposed on 
the upper boundary. At the lower boundary, the vertical component of the velocity is 
set to zero and horizontal components of the turbulent stresses are defined based on the 
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mean logarithmic wind profile assumption as follows (Moeng 1984) 

Ti3 = ~ (^J£ 1 ( 2 - 2 ) 

where z \ , Ui and zq are the vertical spacing, the velocity of the first grid point away from 
the surface, and the roughness height, respectively. This type of boundary condition for 
the surface is a common practice in atmospheric modeling. 

3. Subgrid-scale turbulence modeling 

Turbulent stresses that appear in equation 2.1 are modeled based on the Boussinesq 
eddy viscosity assumption. The Smagorinsky eddy viscosity model (Smagorinsky 1963) 
is a popular approach to represent turbulent stresses 

Tij = -2Cl 2 \S\S ijf (3.1) 

where \S\ is the magnitude of the filtered strain rate tensor, it is defined as 

|£| = yft Sij = - + g£j ■ ( 3 - 2 ) 

In LES, the filter width A = ( dx ■ dy ■ dz) 1 ^ 3 is chosen as the length scale l in the eddy 
viscosity definition (equation 3.1). In the original Smagorinsky model, the dimensionless 
parameter C is an empirical constant, whereas in the dynamic Smagorinsky model (Ger- 
mano et al. 1991), this parameter is computed dynamically during the solution, making 
it a function of space and time. 

In the vicinity of the surface, the numerical resolution is generally limited and near- 
surface models are adopted to model the turbulent stresses properly. We consider three 
near-surface models, and provide a brief explanation of them in the following. 

3.1. Mason and Thomson Model 

This is an arbitrary and simple modification of the length scale in the eddy viscosity 
definition given by equation 3.1. Mason & Thomson (1992) have suggested using the 
following form as the length scale. 


P A 2 + {kz) 2 ^ 

where k is the von Karman constant with a value of 0.41, and 2 is the distance from the 
surface. 


3.2. Canopy Stress Model 

Canopy models are usually adopted to model the structure of the flow within the canopy. 
Brown et al. (2001) have suggested using them as near-surface models in the case of coarse 
numerical resolution. In particular, we follow the implementation of Kirkpatrick et al. 
(2003). Canopy stresses near the surface are formulated as 


f h 

Ti 3 = - / C s COS 

Jo 



\u\uidz, 


(3.4) 


where C s is a normalizing parameter. Specifically, C s is the ratio of the surface stress at 
the bottom to the canopy stress at the first grid point. The canopy height h is taken to 
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be twice the streamwise grid spacing (2 dx) . Canopy stresses are added onto the turbulent 
stresses that are modeled with the dynamic Smagorinsky model. 

3.3. Hybrid RANS/LES model 

We adopt a hybrid RANS/LES approach as the third model. The motivation behind 
this model is that the surface stress boundary conditions represent ensemble-averaged 
quantities, because they are derived from the mean logarithmic wind profile. Hence, it is 
more meaningful to adopt a RANS type eddy viscosity formulation near the surface. 

One can analytically show that the logarithmic wind profile can be derived from the 
following eddy viscosity definition (Panofsky & Dutton 1984) 


vt = ku*z, (3.5) 

where u* is the friction velocity. The above form is also identical to the following definition 

Vt = ( KZ ) 2 ^ = ( kz ) 2 \S\- (3-6) 

Here, U is magnitude of the velocity vector. We apply the above definition in the vicinity 
of the surface, which guarantees a logarithmic wind profile, and switch to the dynamic 
Smagorinsky model away from the surface with the help of a blending function. This 
hybrid RANS/LES model can be written as 

vt = [(1 - exp (-z/h)) 2 ■ (CA 2 ) -I- exp(-z/h) 2 (Kz) 2 ]\S\ (3.7) 

Following the suggestion of van Driest (1956) an exponential form is selected as the blend- 
ing function. Here h is the altitude, where the logarithmic wind profile is not expected 
further. As a rule of thump, a logarithmic velocity profile is expected within the surface 
layer, which is the bottom 10 percent of the atmospheric boundary layer height (100 m 
-200 m) (Stull 1988). h = 150m is used for the present computations. 

4. Numerical schemes 

We use the LES atmospheric research code (DHARMA) for the present simulations. 
The numerical method adopted in DHARMA is described in detail in Stevens & Brether- 
ton (1996) and Stevens et al. (2000). The governing equations are integrated using 
a forward-in-time projection method based on an explicit second-order Runge-Kutta 
scheme (Bell & Marcus 1992). The spatial discretization is performed on a staggered 
grid. A third-order accurate upwind-biased monotonic scheme is used for the advection 
terms, whereas diffusion and pressure gradient terms are discretized using second-order 
accurate central differencing schemes. A direct solver (FFT) is utilized for solving the 
pressure Poisson equation, and the code has been parallelized to run on various platforms 
using MPI. 

5. Results 

In this section, we present results from the simulation of neutrally stratified atmo- 
spheric boundary layer. The computational domain size is 3000 m x 1500 m x 1500 
m with 64 x 32 x 64 grid points in x, y, and z directions, respectively. In the vertical 
direction, the grid points are clustered exponentially near the surface, where the first 
grid point for the horizontal component of the velocity is located at 4 m from the sur- 
face. The flow is driven by the Coriolis force that would balance a 10 ms -1 geostrophic 
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Figure 1 . Comparison of mean wind profile. : log- law, — ¥ — : no near-surface model, 

— o — : hybrid RANS/LES model, — o — : canopy stress model, — x — : Mason-Thomson 
model 


wind in the x direction. Coriolis parameter is set equal to 10 -4 s _1 , and the roughness 
parameter zg has a value of 0.1m. A dimensionless time unit can be defined based on the 
Coriolis parameter(l//). The simulations were run over a period of 10/ _1 , and statistics 
are collected during the time period of last the 4/ _1 . Ensemble averaged vertical profiles 
have been obtained by collecting data at every 6 x 10 -3 dimensionless units (60secs), 
and averaging them both in time and in horizontal space. 

Andren et al. (1994) compared the performance of different LES computer codes to 
simulate a neutrally stratified atmospheric boundary layer. They concluded that results 
were more sensitive to the SGS model in the lower third of the boundary layer, and 
commonly used SGS models have failed to reproduce the logarithmic wind profile in the 
vicinity of the surface. Hence, satisfying the logarithmic wind profile near the surface has 
become one of the measures of good performance for an SGS model. 

Figure 1 shows a comparison of the mean wind profiles for all the models considered 
in the present study. When no near-surface model is used, the dynamic Smagorinsky 
model produces an erroneous profile. This indicates that providing only the surface stress 
boundary condition is not sufficient to model the SGS motions near the surface in the 
case of high Reynolds number and coarse numerical resolution. This observation is in 
agreement with the findings of Cabot et al. (1999). It should be noted that the erroneous 
profile is not surprising because the essence of the dynamic Smagorinsky model is based 
on resolving the inertial subrange. This has not been satisfied with the current coarse 
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Figure 2. Comparison of streamwise component of the total (left) and resolved (right) stress. 

• • • : no near-surface model, : hybrid RANS/LES model, : canopy stress model, 

— — : Mason-Thomson model 

resolution. Hence, SGS models need augmentation to represent the unresolved scales near 
the surface. 

The near-surface models considered in the present study, namely, Mason-Thomson 
model, the canopy stress model and the hybrid RANS/LES model are formulated to 
parameterize the unresolved scales near the surface. Both the canopy stress model and the 
hybrid RANS/LES model adopt the dynamic Smagorinsky model away from the surface, 
whereas the Mason-Thomson model adopts the original Smagorinsky model. As seen from 
figure 1, results obtained from these three near-surface models are much closer to the log- 
law profile. However, considerable differences exist among them. The best agreement with 
the log-law near the surface is achieved with the hybrid RANS/LES model, whereas the 
Mason and Thomson model gives the worst agreement. The hybrid RANS/LES model has 
the closest agreement, because the length scale is taken to be the distance from the surface 
in the model, which is known to give the log-law analytically. The canopy stress model 
follows the log-law close to the surface, but a significant deviation is observed further 
from the surface. This indicates that, canopy stresses should be distributed over a larger 
distance from the surface. We reiterate that, based on field experiments, a logarithmic 
velocity profile is expected in the surface layer, which is approximately the bottom 10 
percent of the atmospheric boundary layer height (100 m -200 m). 

The streamwise (( uw/ul ) ) component of the total and resolved turbulent stress profiles 
are compared in figure 2. By total stress, we mean the resolved turbulent stresses plus the 
modeled SGS turbulent stresses. As seen from this plot, simulations with surface models 
predict an almost linear variation of the total turbulent stress in agreement with surface 
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Figure 3. Comparison of spanwise component of the total (left) and resolved (right) stress. • • • : 

no near-surface model, : hybrid RANS/LES model, : canopy stress model, : 

Mason-Thomson model 


similarity theory (Garratt 1992). However, without a near-surface model, the dynamic 
Smagorinsky model does not predict a linear variation of the total turbulent stress and 
produces wiggles at the surface, indicating the need for near-surface modeling. 

The resolved turbulent stresses are also shown in figure 2. Within the lower one fourth 
of the boundary layer, a large range of turbulent scales are modeled and not resolved 
with the hybrid RANS/LES model, whereas other surface models are effective only within 
the lower one tenth of the boundary layer. The extent of RANS modeling in the hybrid 
RANS/LES model is adjustable through the value of h in equation 3.7. We have found 
that h = 150m is reasonable in terms of satisfying the logarithmic wind profile in accor- 
dance with experimental observations. 

The differences among the three surface models are more pronounced in the spanwise 
({v'w'/v%)) component of the total stress as shown in 3. This spanwise component of the 
total stress is generated due to the rotation of the Earth. The variation of the velocity 
vector with height is known as the Ekman spiral (Stull 1988). 

As can be seen from figure 4, the spanwise component of the velocity is also affected 
because of differences in the distribution of the spanwise component of the total turbulent 
stress. With all three surface models, we see that the spanwise component of the velocity 
does not vanish at the top, whereas it vanishes in the case of the dynamic Smagorinsky 
model without a near-surface model. Apparently, nonvanishing spanwise velocity at the 
top is an outcome of employing the surface models. Since ensemble averaged quantities 
are obtained near the wall, fluctuations of the velocity field are damped due to the action 
of turbulent eddy viscosity or canopy stresses, effectively resulting in a laminar flow field. 
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Figure 4. Comparison of spanwise component (v) of the velocity profiles. • • • : no near-surface 

model, : hybrid RANS/LES model, : canopy stress model, : Mason-Thomson 

model 

As a result, this ensemble averaged region acts like an artificial boundary layer for the 
flow field away from the surface that is modeled in the LES sense. Baggett (1998) have 
reported similar issues in LES of channel flow. 

We present turbulent eddy viscosity profiles in figure 5. Away from the surface, iden- 
tical profiles are obtained with the hybrid RANS/LES model and the canopy stress 
model. Near the surface, the eddy viscosity levels are much lower with the canopy stress 
model. This is because the model modifies the stress terms but not the eddy viscosity. 
On the other hand, eddy viscosity level is much higher near the surface with the hybrid 
RANS/LES model. With such high levels of eddy viscosity, turbulent motions near the 
surface are dampened, resulting in an ensemble averaged velocity field. Among all mod- 
els, the Mason-Thomson predicts higher levels of eddy viscosity away from the surface, 
because it adopts the original Smagorinsky model in that region. 


6. Summary and conclusions 

In this study, we have performed large-eddy simulations of a neutrally stratified at- 
mospheric boundary layer. Particularly, we have assessed the performance of different 
near-surface models in predicting the mean flow structure. The surface of our atmo- 
sphere is always rough and typical numerical resolutions that are utilized in atmospheric 
boundary layer modeling are too coarse near the surface to perform well-resolved large- 
eddy simulations. This necessitates additional modeling in the vicinity of the surface to 
take into account the unresolved turbulent flow, covering a large range of scales. We have 
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Figure 5. Comparison of turbulent eddy viscosity profiles. ■ • ■ : no near-surface model, : 

hybrid RANS/LES model, : canopy stress model, : Mason-Thomson model 

considered the near-surface models of Mason & Thomson (1992), Brown et al. (2001) and 
also suggested a hybrid RANS/LES model that blends Prandtl’s mixing length model 
with the dynamic Smagorinsky model. We have also considered the dynamic Smagorin- 
sky model with no near-surface model to highlight the need for surface modeling in the 
case of coarse numerical resolution. 

Off all the near-surface models considered, the hybrid RANS/LES model shows better 
agreement with mean logarithmic wind profile. In the canopy stress model, it appears that 
surface stresses should be distributed to a larger distance from the surface to improve the 
predictions. The downside of adopting any of the near-surface models is that the spanwise 
component of the velocity vector does not vanish at the top of the computational domain. 

For the current numerical resolution with the hybrid RANS/LES model, the contribution 
from RANS modeling is effective on about the lower one fourth of the domain. A finer 
numerical resolution might help confine the RANS contribution very close to the surface 
and get better predictions at the same time. 

Although comparable predictions can be obtained by either adopting the hybrid RANS/LES 
or the canopy stress model, the hybrid RANS/LES model has the advantage of easy im- 
plementation. In a parallel effort (Senocak et al. 2004), we are implementing the immersed 
boundary method in DHARMA. We find that it is much more convenient to implement 
the hybrid RANS/LES model to simulate flow over topography with the immersed bound- 
ary method because only the eddy viscosity term needs modification. However, it might 
not be straightforward to adopt the canopy stress model in Cartesian grid techniques for 
complex topography simulations. 
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